"""
DaoMath Condensed Matter Physics (道数凝聚态物理)
从 ±0 先天均衡构造凝聚态理论

凝聚态物理的 ±0 诠释:
    电子 ↔ 空穴 (对偶)
    导体 ↔ 绝缘体 (对偶)
    正常态 ↔ 超导态 (对偶)
    长程序 ↔ 无序 (对偶)
    
理论层级:
    ±0 (先天均衡)
     ↓
    晶格结构 (Lattice)
     ↓
    能带理论 (Band Structure)
     ↓
    准粒子激发
     ↓
    多体效应
     ↓
    相变与对称破缺
     ↓
    超导与超流
     ↓
    拓扑相
"""

from __future__ import annotations
from typing import List, Tuple, Callable, Optional, Dict, Any, Union
from dataclasses import dataclass, field
from abc import ABC, abstractmethod
from enum import Enum
import numpy as np
import scipy.linalg as la
from scipy.integrate import quad, odeint
from scipy.optimize import minimize, fsolve
from scipy.special import ellipk, ellipe
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

from core.axioms import InnateBalance, DualPair


# ============================================================================
# 物理常数
# ============================================================================

class CMConstants:
    """凝聚态物理常数"""
    HBAR = 1.054571817e-34      # ℏ (J·s)
    M_E = 9.1093837015e-31      # 电子质量 (kg)
    E_CHARGE = 1.602176634e-19  # 基本电荷 (C)
    K_B = 1.380649e-23          # Boltzmann 常数 (J/K)
    MU_B = 9.2740100783e-24     # Bohr 磁子 (J/T)
    
    # 简化单位制
    HBAR_AU = 1.0
    M_E_AU = 1.0
    E_CHARGE_AU = 1.0
    K_B_AU = 1.0


# ============================================================================
# Level 1: 晶格结构
# ============================================================================

class Bravais:
    """
    Bravais 晶格类型
    """
    SIMPLE_CUBIC = "sc"
    BODY_CENTERED_CUBIC = "bcc"
    FACE_CENTERED_CUBIC = "fcc"
    HEXAGONAL = "hex"
    SQUARE = "square"
    HONEYCOMB = "honeycomb"


class LatticeVector:
    """
    晶格矢量
    
    基矢: a₁, a₂, a₃
    """
    
    def __init__(self, vectors: List[np.ndarray]):
        """
        Args:
            vectors: 基矢列表 [a₁, a₂, a₃] 或 [a₁, a₂]
        """
        self.vectors = [np.array(v) for v in vectors]
        self.dim = len(vectors)
        
        # 验证线性独立
        if self.dim > 1:
            matrix = np.column_stack(self.vectors)
            if abs(np.linalg.det(matrix[:self.dim, :self.dim])) < 1e-10:
                raise ValueError("Lattice vectors must be linearly independent")
    
    def primitive_cell_volume(self) -> float:
        """
        原胞体积
        
        V = |a₁ · (a₂ × a₃)|  (3D)
        V = |a₁ × a₂|         (2D)
        """
        if self.dim == 1:
            return np.linalg.norm(self.vectors[0])
        elif self.dim == 2:
            return abs(np.cross(self.vectors[0], self.vectors[1]))
        elif self.dim == 3:
            return abs(np.dot(self.vectors[0], 
                            np.cross(self.vectors[1], self.vectors[2])))
        return 0.0
    
    def reciprocal_lattice(self) -> LatticeVector:
        """
        倒格子
        
        b₁ = 2π (a₂ × a₃) / V
        b₂ = 2π (a₃ × a₁) / V
        b₃ = 2π (a₁ × a₂) / V
        
        满足: aᵢ · bⱼ = 2π δᵢⱼ
        """
        V = self.primitive_cell_volume()
        
        if self.dim == 1:
            b1 = 2 * np.pi * self.vectors[0] / V**2
            return LatticeVector([b1])
        
        elif self.dim == 2:
            a1, a2 = self.vectors
            # 2D 倒格子（需要扩展到 3D 计算）
            a1_3d = np.append(a1, 0)
            a2_3d = np.append(a2, 0)
            
            b1 = 2 * np.pi * np.cross(a2_3d, [0, 0, 1]) / V
            b2 = 2 * np.pi * np.cross([0, 0, 1], a1_3d) / V
            
            return LatticeVector([b1[:2], b2[:2]])
        
        elif self.dim == 3:
            a1, a2, a3 = self.vectors
            
            b1 = 2 * np.pi * np.cross(a2, a3) / V
            b2 = 2 * np.pi * np.cross(a3, a1) / V
            b3 = 2 * np.pi * np.cross(a1, a2) / V
            
            return LatticeVector([b1, b2, b3])
        
        raise ValueError("Invalid dimension")
    
    @staticmethod
    def simple_cubic(a: float) -> LatticeVector:
        """简单立方晶格"""
        return LatticeVector([
            [a, 0, 0],
            [0, a, 0],
            [0, 0, a]
        ])
    
    @staticmethod
    def fcc(a: float) -> LatticeVector:
        """面心立方晶格"""
        return LatticeVector([
            [0, a/2, a/2],
            [a/2, 0, a/2],
            [a/2, a/2, 0]
        ])
    
    @staticmethod
    def bcc(a: float) -> LatticeVector:
        """体心立方晶格"""
        return LatticeVector([
            [-a/2, a/2, a/2],
            [a/2, -a/2, a/2],
            [a/2, a/2, -a/2]
        ])
    
    @staticmethod
    def square(a: float) -> LatticeVector:
        """正方晶格 (2D)"""
        return LatticeVector([
            [a, 0],
            [0, a]
        ])
    
    @staticmethod
    def honeycomb(a: float) -> LatticeVector:
        """
        蜂窝晶格 (石墨烯)
        
        两个三角子格的组合
        """
        return LatticeVector([
            [a, 0],
            [a/2, a*np.sqrt(3)/2]
        ])


class CrystalStructure:
    """
    晶体结构
    
    = 晶格 + 基
    
    基于 ±0 的理解:
    - 周期性：平移对称（均衡）
    - 基：不同原子类型的叠加
    """
    
    def __init__(self, 
                 lattice: LatticeVector,
                 basis: List[np.ndarray],
                 atomic_species: Optional[List[str]] = None):
        """
        Args:
            lattice: 晶格矢量
            basis: 基位置（相对原胞）
            atomic_species: 原子种类
        """
        self.lattice = lattice
        self.basis = [np.array(b) for b in basis]
        self.num_basis = len(basis)
        self.species = atomic_species or ['X'] * self.num_basis
    
    def atom_position(self, 
                     lattice_index: Tuple[int, ...],
                     basis_index: int) -> np.ndarray:
        """
        原子位置
        
        R = n₁a₁ + n₂a₂ + n₃a₃ + r_basis
        """
        position = np.zeros(self.lattice.dim)
        
        for i, n in enumerate(lattice_index):
            position += n * self.lattice.vectors[i]
        
        position += self.basis[basis_index]
        
        return position
    
    def nearest_neighbors(self, 
                         max_distance: float = 2.0) -> List[Tuple[np.ndarray, float]]:
        """
        最近邻原子
        
        Returns:
            [(相对位置, 距离)]
        """
        neighbors = []
        
        # 搜索范围
        search_range = range(-2, 3)
        
        for n1 in search_range:
            for n2 in search_range:
                for n3 in search_range if self.lattice.dim == 3 else [0]:
                    
                    lattice_idx = (n1, n2, n3)[:self.lattice.dim]
                    
                    for basis_idx in range(self.num_basis):
                        pos = self.atom_position(lattice_idx, basis_idx)
                        distance = np.linalg.norm(pos)
                        
                        if 0 < distance < max_distance:
                            neighbors.append((pos, distance))
        
        # 按距离排序
        neighbors.sort(key=lambda x: x[1])
        
        return neighbors
    
    def coordination_number(self, tolerance: float = 0.1) -> int:
        """
        配位数（最近邻数目）
        """
        neighbors = self.nearest_neighbors()
        
        if not neighbors:
            return 0
        
        min_distance = neighbors[0][1]
        
        count = sum(1 for _, d in neighbors 
                   if abs(d - min_distance) < tolerance)
        
        return count
    
    @staticmethod
    def diamond(a: float) -> CrystalStructure:
        """
        金刚石结构
        
        FCC + 基: (0,0,0) 和 (a/4, a/4, a/4)
        """
        lattice = LatticeVector.fcc(a)
        basis = [
            [0, 0, 0],
            [a/4, a/4, a/4]
        ]
        
        return CrystalStructure(lattice, basis, ['C', 'C'])
    
    @staticmethod
    def graphene(a: float) -> CrystalStructure:
        """
        石墨烯
        
        蜂窝晶格 + 基: A 和 B 原子
        """
        lattice = LatticeVector.honeycomb(a)
        
        # 两个原子位置
        delta = a / np.sqrt(3)
        basis = [
            [0, 0],
            [delta/2, delta*np.sqrt(3)/2]
        ]
        
        return CrystalStructure(lattice, basis, ['A', 'B'])


# ============================================================================
# Level 2: 能带理论
# ============================================================================

class BrillouinZone:
    """
    Brillouin 区
    
    倒空间中的 Wigner-Seitz 原胞
    """
    
    def __init__(self, reciprocal_lattice: LatticeVector):
        """
        Args:
            reciprocal_lattice: 倒格子
        """
        self.reciprocal = reciprocal_lattice
        self.dim = reciprocal_lattice.dim
    
    def high_symmetry_points(self) -> Dict[str, np.ndarray]:
        """
        高对称点
        
        例: Γ (0,0,0), X, M, K 等
        """
        if self.dim == 2:
            # 正方晶格
            return {
                'Γ': np.array([0, 0]),
                'X': np.array([np.pi, 0]),
                'M': np.array([np.pi, np.pi])
            }
        elif self.dim == 3:
            # 简单立方
            return {
                'Γ': np.array([0, 0, 0]),
                'X': np.array([np.pi, 0, 0]),
                'M': np.array([np.pi, np.pi, 0]),
                'R': np.array([np.pi, np.pi, np.pi])
            }
        
        return {'Γ': np.zeros(self.dim)}
    
    def k_path(self, 
              points: List[str],
              num_points: int = 100) -> Tuple[np.ndarray, List[str]]:
        """
        高对称路径
        
        Args:
            points: 高对称点标签列表，如 ['Γ', 'X', 'M', 'Γ']
        
        Returns:
            (k 点数组, 标签位置)
        """
        sym_points = self.high_symmetry_points()
        
        k_points = []
        labels = []
        
        for i in range(len(points) - 1):
            start = sym_points[points[i]]
            end = sym_points[points[i+1]]
            
            # 线性插值
            segment = np.linspace(start, end, num_points)
            k_points.append(segment)
            
            if i == 0:
                labels.append((0, points[i]))
            labels.append((len(k_points[-1]) - 1, points[i+1]))
        
        k_path_array = np.vstack(k_points)
        
        return (k_path_array, labels)


class TightBindingModel:
    """
    紧束缚模型 (Tight-Binding)
    
    H = Σᵢ εᵢ |i⟩⟨i| + Σ_{<i,j>} tᵢⱼ |i⟩⟨j|
    
    基于 ±0 的理解:
    - 在位能 εᵢ: 局域能量
    - 跃迁积分 tᵢⱼ: 电子在格点间的跃迁（正负关联）
    """
    
    def __init__(self,
                 crystal: CrystalStructure,
                 hopping: Dict[Tuple[int, int], float],
                 onsite_energy: float = 0.0):
        """
        Args:
            crystal: 晶体结构
            hopping: 跃迁积分 {(i, j): t_ij}
            onsite_energy: 在位能
        """
        self.crystal = crystal
        self.hopping = hopping
        self.epsilon = onsite_energy
    
    def hamiltonian(self, k: np.ndarray) -> np.ndarray:
        """
        Bloch Hamiltonian
        
        H(k) = Σᵢⱼ h_ij e^{ik·(Rᵢ-Rⱼ)}
        
        Args:
            k: 晶体动量
        
        Returns:
            Hamiltonian 矩阵
        """
        n_basis = self.crystal.num_basis
        H = np.zeros((n_basis, n_basis), dtype=complex)
        
        # 在位能
        for i in range(n_basis):
            H[i, i] = self.epsilon
        
        # 跃迁项
        neighbors = self.crystal.nearest_neighbors()
        
        for (i, j), t_ij in self.hopping.items():
            # 计算相位因子
            for R, _ in neighbors:
                phase = np.exp(1j * np.dot(k, R))
                H[i, j] += t_ij * phase
        
        return H
    
    def band_structure(self, 
                      k_path: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
        """
        能带结构
        
        E_n(k): 第 n 条能带
        
        Returns:
            (k 点, 本征能量)
        """
        num_k = len(k_path)
        num_bands = self.crystal.num_basis
        
        energies = np.zeros((num_k, num_bands))
        
        for i, k in enumerate(k_path):
            H_k = self.hamiltonian(k)
            eigenvalues = np.linalg.eigvalsh(H_k)
            energies[i] = eigenvalues
        
        return (k_path, energies)
    
    def density_of_states(self,
                         energy_range: Tuple[float, float],
                         num_k: int = 50,
                         num_energy: int = 100) -> Tuple[np.ndarray, np.ndarray]:
        """
        态密度 (DOS)
        
        ρ(E) = Σ_nk δ(E - E_n(k))
        
        Returns:
            (能量, 态密度)
        """
        E_min, E_max = energy_range
        energies = np.linspace(E_min, E_max, num_energy)
        dos = np.zeros(num_energy)
        
        # 生成 k 点网格
        if self.crystal.lattice.dim == 1:
            k_grid = np.linspace(-np.pi, np.pi, num_k)
        elif self.crystal.lattice.dim == 2:
            k_grid = np.array([[kx, ky] 
                              for kx in np.linspace(-np.pi, np.pi, num_k)
                              for ky in np.linspace(-np.pi, np.pi, num_k)])
        else:
            k_grid = np.array([[kx, ky, kz]
                              for kx in np.linspace(-np.pi, np.pi, num_k)
                              for ky in np.linspace(-np.pi, np.pi, num_k)
                              for kz in np.linspace(-np.pi, np.pi, num_k)])
        
        # 计算所有 k 点的能量
        bandwidth = 0.1  # Gaussian 展宽
        
        for k in k_grid:
            H_k = self.hamiltonian(k)
            E_k = np.linalg.eigvalsh(H_k)
            
            # 贡献到 DOS
            for E in E_k:
                dos += np.exp(-(energies - E)**2 / (2 * bandwidth**2))
        
        # 归一化
        dos /= (np.sqrt(2 * np.pi) * bandwidth * len(k_grid))
        
        return (energies, dos)
    
    def fermi_surface(self, 
                     E_F: float,
                     num_k: int = 50) -> List[np.ndarray]:
        """
        Fermi 面
        
        能带与 Fermi 能级的交线
        
        Returns:
            Fermi 面上的 k 点列表
        """
        if self.crystal.lattice.dim != 2:
            raise NotImplementedError("Only 2D Fermi surface implemented")
        
        fermi_points = []
        
        k_range = np.linspace(-np.pi, np.pi, num_k)
        
        for kx in k_range:
            for ky in k_range:
                k = np.array([kx, ky])
                H_k = self.hamiltonian(k)
                eigenvalues = np.linalg.eigvalsh(H_k)
                
                # 检查是否穿过 Fermi 能级
                if np.any(np.abs(eigenvalues - E_F) < 0.1):
                    fermi_points.append(k)
        
        return fermi_points
    
    @staticmethod
    def graphene_model(t: float = 1.0) -> TightBindingModel:
        """
        石墨烯紧束缚模型
        
        最近邻跃迁
        """
        graphene = CrystalStructure.graphene(1.0)
        
        # A-B 原子间跃迁
        hopping = {
            (0, 1): -t,  # A → B
            (1, 0): -t   # B → A
        }
        
        return TightBindingModel(graphene, hopping)


class BandGap:
    """
    能隙
    
    导带底与价带顶的能量差
    
    基于 ±0 的理解:
    - 直接带隙: 同一 k 点（对称）
    - 间接带隙: 不同 k 点（破缺）
    """
    
    def __init__(self, 
                 E_g: float,
                 gap_type: str = 'direct',
                 k_valence: Optional[np.ndarray] = None,
                 k_conduction: Optional[np.ndarray] = None):
        """
        Args:
            E_g: 能隙大小
            gap_type: 'direct' 或 'indirect'
            k_valence: 价带顶 k 点
            k_conduction: 导带底 k 点
        """
        self.E_g = E_g
        self.gap_type = gap_type
        self.k_v = k_valence
        self.k_c = k_conduction
    
    def is_insulator(self, temperature: float = 0.0) -> bool:
        """判断是否为绝缘体"""
        k_B = CMConstants.K_B_AU
        
        # 室温下能隙 > 3 eV 为绝缘体
        # 否则为半导体
        if self.E_g > 3.0:
            return True
        
        # 考虑热激发
        if temperature > 0:
            thermal_energy = k_B * temperature
            return self.E_g > 10 * thermal_energy
        
        return self.E_g > 0
    
    def optical_absorption_threshold(self) -> float:
        """
        光吸收阈值能量
        
        直接带隙: E_threshold = E_g
        间接带隙: 需要声子辅助
        """
        if self.gap_type == 'direct':
            return self.E_g
        else:
            # 间接跃迁需要考虑声子能量
            return self.E_g  # 简化


class EffectiveMass:
    """
    有效质量
    
    m* = ℏ² / (∂²E/∂k²)
    
    体现能带曲率
    """
    
    @staticmethod
    def calculate(band_energy: Callable[[np.ndarray], float],
                 k_point: np.ndarray,
                 direction: int = 0,
                 hbar: float = CMConstants.HBAR_AU) -> float:
        """
        计算有效质量
        
        Args:
            band_energy: E(k) 函数
            k_point: 计算点
            direction: k 空间方向
        """
        dk = 1e-4
        
        k_plus = k_point.copy()
        k_plus[direction] += dk
        
        k_minus = k_point.copy()
        k_minus[direction] -= dk
        
        # 二阶导数
        E_plus = band_energy(k_plus)
        E_center = band_energy(k_point)
        E_minus = band_energy(k_minus)
        
        d2E_dk2 = (E_plus - 2*E_center + E_minus) / dk**2
        
        if abs(d2E_dk2) < 1e-10:
            return float('inf')  # 平带
        
        return hbar**2 / d2E_dk2


# ============================================================================
# Level 3: 电子-声子相互作用
# ============================================================================

class Phonon:
    """
    声子（晶格振动量子）
    
    准粒子激发
    
    基于 ±0 的理解:
    - 声学支：长波振动
    - 光学支：短波振动
    - 正负位移的对偶
    """
    
    def __init__(self,
                 crystal: CrystalStructure,
                 force_constants: np.ndarray):
        """
        Args:
            crystal: 晶体结构
            force_constants: 力常数矩阵
        """
        self.crystal = crystal
        self.K = force_constants
    
    def dynamical_matrix(self, k: np.ndarray) -> np.ndarray:
        """
        动力学矩阵
        
        D(k) = Σ_R K(R) e^{ik·R} / √(M_i M_j)
        """
        # 简化实现
        n = self.crystal.num_basis
        D = np.zeros((3*n, 3*n), dtype=complex)
        
        # Fourier 变换力常数
        for i in range(3*n):
            for j in range(3*n):
                D[i, j] = self.K[i, j]  # 简化，不考虑相位
        
        return D
    
    def phonon_dispersion(self, 
                         k_path: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
        """
        声子色散关系
        
        ω_λ(k): 第 λ 支声子频率
        
        Returns:
            (k 点, 频率)
        """
        num_k = len(k_path)
        num_modes = 3 * self.crystal.num_basis
        
        frequencies = np.zeros((num_k, num_modes))
        
        for i, k in enumerate(k_path):
            D_k = self.dynamical_matrix(k)
            eigenvalues = np.linalg.eigvalsh(D_k)
            
            # ω² = eigenvalue, 取正根
            frequencies[i] = np.sqrt(np.abs(eigenvalues))
        
        return (k_path, frequencies)
    
    def debye_temperature(self, 
                         sound_velocity: float,
                         density: float) -> float:
        """
        Debye 温度
        
        θ_D = ℏω_D / k_B
        
        其中 ω_D = v_s (6π²n)^(1/3)
        """
        k_B = CMConstants.K_B_AU
        hbar = CMConstants.HBAR_AU
        
        omega_D = sound_velocity * (6 * np.pi**2 * density)**(1/3)
        
        return hbar * omega_D / k_B


class ElectronPhononCoupling:
    """
    电子-声子耦合
    
    H_ep = Σ_kqλ g_λ(k,q) c†_{k+q} c_k (b_qλ + b†_{-qλ})
    
    导致：
    - 电阻率（散射）
    - 超导（Cooper 对）
    """
    
    def __init__(self, coupling_constant: float):
        """
        Args:
            coupling_constant: 耦合强度 λ
        """
        self.lambda_ep = coupling_constant
    
    def electron_lifetime(self,
                         temperature: float,
                         debye_temp: float) -> float:
        """
        电子寿命（由声子散射决定）
        
        1/τ ∝ T (T < θ_D)
        1/τ ∝ const (T > θ_D)
        """
        if temperature < debye_temp:
            return 1.0 / (self.lambda_ep * temperature)
        else:
            return 1.0 / (self.lambda_ep * debye_temp)
    
    def resistivity(self, temperature: float) -> float:
        """
        电阻率
        
        ρ(T) ∝ 1/τ
        """
        tau = self.electron_lifetime(temperature, 300)  # 假设 θ_D = 300K
        
        return 1.0 / tau


# ============================================================================
# Level 4: 超导理论
# ============================================================================

class BCSTheory:
    """
    BCS 超导理论
    
    Cooper 对：两个电子通过声子交换形成配对
    
    基于 ±0 的理解:
    - 配对态: (k↑, -k↓) 动量相反（对偶）
    - 能隙: 超导态与正常态的能量差
    - 对称破缺: U(1) 规范对称自发破缺
    """
    
    def __init__(self,
                 coupling_strength: float,
                 debye_frequency: float,
                 fermi_energy: float):
        """
        Args:
            coupling_strength: BCS 耦合 λ = N(0)V
            debye_frequency: Debye 频率 ω_D
            fermi_energy: Fermi 能量 E_F
        """
        self.lambda_bcs = coupling_strength
        self.omega_D = debye_frequency
        self.E_F = fermi_energy
    
    def critical_temperature(self) -> float:
        """
        临界温度
        
        T_c = 1.14 θ_D exp(-1/λ)
        
        弱耦合 BCS 公式
        """
        k_B = CMConstants.K_B_AU
        hbar = CMConstants.HBAR_AU
        
        theta_D = hbar * self.omega_D / k_B
        
        if self.lambda_bcs <= 0:
            return 0.0
        
        T_c = 1.14 * theta_D * np.exp(-1.0 / self.lambda_bcs)
        
        return T_c
    
    def energy_gap(self, temperature: float) -> float:
        """
        超导能隙
        
        Δ(T): 温度依赖的能隙
        
        Δ(0) = 1.76 k_B T_c (T=0)
        """
        k_B = CMConstants.K_B_AU
        T_c = self.critical_temperature()
        
        if temperature >= T_c:
            return 0.0  # 正常态
        
        # BCS 温度依赖（简化）
        Delta_0 = 1.76 * k_B * T_c
        
        # 近似公式
        Delta_T = Delta_0 * np.sqrt(1 - (temperature / T_c)**2)
        
        return Delta_T
    
    def coherence_length(self, temperature: float = 0.0) -> float:
        """
        相干长度
        
        ξ_0 = ℏv_F / (π Δ(0))
        
        Cooper 对的空间尺度
        """
        hbar = CMConstants.HBAR_AU
        
        # Fermi 速度
        m_e = CMConstants.M_E_AU
        v_F = np.sqrt(2 * self.E_F / m_e)
        
        Delta = self.energy_gap(temperature)
        
        if Delta == 0:
            return float('inf')
        
        return hbar * v_F / (np.pi * Delta)
    
    def penetration_depth(self, temperature: float = 0.0) -> float:
        """
        穿透深度（London 穿透深度）
        
        λ_L = √(m / (μ_0 n_s e²))
        
        磁场穿入超导体的深度
        """
        m_e = CMConstants.M_E_AU
        e = CMConstants.E_CHARGE_AU
        
        # 超导电子密度（简化）
        T_c = self.critical_temperature()
        if temperature >= T_c:
            return float('inf')
        
        # 温度依赖（简化）
        n_s_0 = 1e28  # 参考密度
        n_s = n_s_0 * (1 - (temperature / T_c)**4)
        
        mu_0 = 4 * np.pi * 1e-7  # 真空磁导率
        
        lambda_L = np.sqrt(m_e / (mu_0 * n_s * e**2))
        
        return lambda_L
    
    def critical_field(self, temperature: float) -> float:
        """
        临界磁场
        
        H_c(T) = H_c(0) [1 - (T/T_c)²]
        
        超过此磁场，超导态被破坏
        """
        T_c = self.critical_temperature()
        
        if temperature >= T_c:
            return 0.0
        
        # H_c(0) = Φ_0 / (2π λ_L ξ)
        Phi_0 = 2.067833848e-15  # 磁通量子 (Wb)
        
        lambda_L = self.penetration_depth(0)
        xi = self.coherence_length(0)
        
        H_c_0 = Phi_0 / (2 * np.pi * lambda_L * xi)
        
        H_c = H_c_0 * (1 - (temperature / T_c)**2)
        
        return H_c
    
    def dos_superconductor(self, 
                          energy: np.ndarray,
                          temperature: float = 0.0) -> np.ndarray:
        """
        超导态密度
        
        N_s(E) = N_n(E) |E| / √(E² - Δ²)  (|E| > Δ)
        N_s(E) = 0                         (|E| < Δ)
        
        体现 ±0: 能隙内态密度为零（禁区）
        """
        Delta = self.energy_gap(temperature)
        
        # 正常态密度（常数）
        N_n = 1.0
        
        dos = np.zeros_like(energy)
        
        for i, E in enumerate(energy):
            if abs(E) > Delta:
                dos[i] = N_n * abs(E) / np.sqrt(E**2 - Delta**2)
        
        return dos


class GinzburgLandau:
    """
    Ginzburg-Landau 理论
    
    唯象理论，描述超导相变
    
    自由能: F = ∫ dx [α|ψ|² + β|ψ|⁴/2 + |(-iℏ∇ - eA)ψ|²/(2m)]
    
    基于 ±0 的理解:
    - 序参量 ψ: 超导态幅度
    - α < 0: 超导相（自发对称破缺）
    - α > 0: 正常相
    """
    
    def __init__(self, 
                 alpha: float,
                 beta: float,
                 coherence_length: float):
        """
        Args:
            alpha: 线性系数（温度依赖）
            beta: 非线性系数
            coherence_length: 相干长度 ξ
        """
        self.alpha = alpha
        self.beta = beta
        self.xi = coherence_length
    
    def order_parameter(self) -> complex:
        """
        平衡序参量
        
        |ψ_0|² = -α/β  (α < 0)
        ψ_0 = 0        (α > 0)
        """
        if self.alpha >= 0:
            return 0.0 + 0.0j
        
        psi_magnitude = np.sqrt(-self.alpha / self.beta)
        
        # 任意相位（自发对称破缺）
        phase = 0.0
        
        return psi_magnitude * np.exp(1j * phase)
    
    def free_energy_density(self, psi: complex) -> float:
        """
        自由能密度
        
        f = α|ψ|² + β|ψ|⁴/2
        
        （不含梯度项）
        """
        psi_squared = abs(psi)**2
        
        return self.alpha * psi_squared + 0.5 * self.beta * psi_squared**2
    
    def condensation_energy(self) -> float:
        """
        凝聚能
        
        ΔF = F_s - F_n = -α²/(2β)
        
        超导态比正常态低的自由能
        """
        if self.alpha >= 0:
            return 0.0
        
        return -self.alpha**2 / (2 * self.beta)


# ============================================================================
# Level 5: 磁性
# ============================================================================

class Magnetism:
    """
    磁性
    
    基于 ±0 的理解:
    - 顺磁: 自旋随机（无序）
    - 铁磁: 自旋平行（有序）
    - 反铁磁: 自旋反平行（交替序）
    """
    pass


class HeisenbergModel:
    """
    Heisenberg 模型
    
    H = -J Σ_{<i,j>} S_i · S_j
    
    描述局域自旋相互作用
    """
    
    def __init__(self, 
                 coupling: float,
                 lattice_size: int,
                 dimension: int = 1):
        """
        Args:
            coupling: 交换积分 J
            lattice_size: 格点数
            dimension: 维度
        """
        self.J = coupling
        self.L = lattice_size
        self.dim = dimension
        
        # 自旋配置（经典自旋）
        self.spins = np.random.randn(self.L, 3)
        
        # 归一化
        for i in range(self.L):
            self.spins[i] /= np.linalg.norm(self.spins[i])
    
    def energy(self) -> float:
        """
        系统能量
        
        E = -J Σ S_i · S_{i+1}
        """
        E = 0.0
        
        for i in range(self.L):
            j = (i + 1) % self.L  # 周期边界
            E -= self.J * np.dot(self.spins[i], self.spins[j])
        
        return E
    
    def magnetization(self) -> np.ndarray:
        """
        总磁矩
        
        M = Σ_i S_i
        """
        return np.sum(self.spins, axis=0)
    
    def is_ferromagnetic(self) -> bool:
        """判断是否为铁磁相"""
        M = self.magnetization()
        M_magnitude = np.linalg.norm(M)
        
        # 磁化强度接近饱和
        return M_magnitude > 0.8 * self.L
    
    def curie_temperature(self, S: float = 0.5) -> float:
        """
        Curie 温度（平均场近似）
        
        T_c = (2d/3) J S(S+1) / k_B
        
        其中 d = 配位数
        """
        k_B = CMConstants.K_B_AU
        
        coordination = 2 * self.dim  # 简化
        
        T_c = (2 * coordination / 3) * self.J * S * (S + 1) / k_B
        
        return T_c


# ============================================================================
# Level 6: 拓扑相
# ============================================================================

class TopologicalPhase:
    """
    拓扑相
    
    基于 ±0 的理解:
    - 拓扑不变量：整数（量子化）
    - 边界态：体态的对偶
    - 拓扑保护：对微扰鲁棒
    """
    pass


class ChernNumber:
    """
    Chern 数（拓扑不变量）
    
    C = (1/2π) ∫∫ F(k) dk_x dk_y
    
    其中 F(k) = ∂_x A_y - ∂_y A_x 是 Berry 曲率
    """
    
    @staticmethod
    def calculate(hamiltonian: Callable[[np.ndarray], np.ndarray],
                 num_k: int = 50) -> int:
        """
        计算 Chern 数
        
        Args:
            hamiltonian: H(k) 函数
            num_k: k 空间采样点数
        
        Returns:
            Chern 数（整数）
        """
        # 创建 k 空间网格
        kx_vals = np.linspace(-np.pi, np.pi, num_k)
        ky_vals = np.linspace(-np.pi, np.pi, num_k)
        
        chern = 0.0
        
        for i in range(num_k - 1):
            for j in range(num_k - 1):
                # 计算一个小方格的贡献
                k_points = [
                    np.array([kx_vals[i], ky_vals[j]]),
                    np.array([kx_vals[i+1], ky_vals[j]]),
                    np.array([kx_vals[i+1], ky_vals[j+1]]),
                    np.array([kx_vals[i], ky_vals[j+1]])
                ]
                
                # Berry 相位
                U_link = 1.0 + 0.0j
                
                for idx in range(4):
                    k1 = k_points[idx]
                    k2 = k_points[(idx + 1) % 4]
                    
                    # 本征态重叠
                    H1 = hamiltonian(k1)
                    H2 = hamiltonian(k2)
                    
                    _, psi1 = np.linalg.eigh(H1)
                    _, psi2 = np.linalg.eigh(H2)
                    
                    # 最低能带
                    overlap = np.vdot(psi1[:, 0], psi2[:, 0])
                    U_link *= overlap / abs(overlap)
                
                # Berry 曲率贡献
                F = -1j * np.log(U_link)
                chern += np.imag(F)
        
        chern /= (2 * np.pi)
        
        return int(round(chern))


class QuantumHallEffect:
    """
    量子霍尔效应
    
    σ_xy = ν e²/h
    
    其中 ν 是填充因子（整数或分数）
    
    体现 ±0: 边界态的手性流动（单向）
    """
    
    def __init__(self, filling_factor: Union[int, float]):
        """
        Args:
            filling_factor: 填充因子 ν
        """
        self.nu = filling_factor
    
    def hall_conductivity(self) -> float:
        """
        霍尔电导
        
        σ_xy = ν e²/h
        """
        e = CMConstants.E_CHARGE
        h = CMConstants.HBAR * 2 * np.pi
        
        return self.nu * e**2 / h
    
    def is_integer(self) -> bool:
        """判断是整数量子霍尔效应还是分数量子霍尔效应"""
        return abs(self.nu - round(self.nu)) < 1e-10


class TopologicalInsulator:
    """
    拓扑绝缘体
    
    体态：绝缘
    边界：导电（受拓扑保护）
    
    Z₂ 不变量
    """
    
    def __init__(self, 
                 bulk_gap: float,
                 surface_states: bool = True):
        """
        Args:
            bulk_gap: 体能隙
            surface_states: 是否有表面态
        """
        self.E_g = bulk_gap
        self.has_surface_states = surface_states
    
    def is_topological(self) -> bool:
        """判断是否为非平凡拓扑绝缘体"""
        return self.E_g > 0 and self.has_surface_states
    
    def z2_invariant(self) -> int:
        """
        Z₂ 拓扑不变量
        
        ν = 0 (平凡) 或 1 (非平凡)
        """
        return 1 if self.is_topological() else 0


# ============================================================================
# 综合演示系统
# ============================================================================

class CondensedMatterDemo:
    """凝聚态物理综合演示"""
    
    @staticmethod
    def demonstrate_lattice():
        """演示晶格结构"""
        print("\n" + "=" * 70)
        print("晶格结构")
        print("=" * 70)
        
        # 简单立方
        sc = LatticeVector.simple_cubic(1.0)
        print(f"\n简单立方:")
        print(f"  原胞体积: {sc.primitive_cell_volume():.4f}")
        
        # 倒格子
        reciprocal = sc.reciprocal_lattice()
        print(f"  倒格子基矢:")
        for i, b in enumerate(reciprocal.vectors):
            print(f"    b_{i+1} = {b}")
        
        # FCC
        fcc = LatticeVector.fcc(1.0)
        print(f"\nFCC:")
        print(f"  原胞体积: {fcc.primitive_cell_volume():.4f}")
        
        # 石墨烯
        graphene = CrystalStructure.graphene(1.42)  # C-C 键长 1.42 Å
        print(f"\n石墨烯:")
        print(f"  基原子数: {graphene.num_basis}")
        print(f"  配位数: {graphene.coordination_number()}")
        
        neighbors = graphene.nearest_neighbors()[:3]
        print(f"  最近邻距离:")
        for pos, dist in neighbors:
            print(f"    d = {dist:.4f}")
    
    @staticmethod
    def demonstrate_band_structure():
        """演示能带结构"""
        print("\n" + "=" * 70)
        print("能带结构")
        print("=" * 70)
        
        # 石墨烯紧束缚模型
        tb = TightBindingModel.graphene_model(t=1.0)
        
        # Brillouin 区
        bz = BrillouinZone(tb.crystal.lattice.reciprocal_lattice())
        
        # 高对称路径
        k_path, labels = bz.k_path(['Γ', 'X', 'M', 'Γ'], num_points=100)
        
        # 能带结构
        k_points, energies = tb.band_structure(k_path)
        
        print(f"\n石墨烯能带:")
        print(f"  能带数: {energies.shape[1]}")
        print(f"  在 Γ 点:")
        E_gamma = energies[0]
        for i, E in enumerate(E_gamma):
            print(f"    E_{i} = {E:.4f}")
        
        # Dirac 点
        print(f"\n  Dirac 点（K 点）:")
        print(f"    线性色散关系")
        
        # 态密度
        dos_E, dos_rho = tb.density_of_states((-3, 3), num_k=20)
        print(f"\n态密度:")
        print(f"  能量范围: [{dos_E[0]:.2f}, {dos_E[-1]:.2f}]")
        print(f"  Fermi 能级处 DOS: {dos_rho[len(dos_rho)//2]:.4f}")
    
    @staticmethod
    def demonstrate_superconductivity():
        """演示超导"""
        print("\n" + "=" * 70)
        print("BCS 超导理论")
        print("=" * 70)
        
        # 典型超导体参数（如铝）
        lambda_bcs = 0.18
        omega_D = 4.2e13  # rad/s
        E_F = 11.7  # eV
        
        bcs = BCSTheory(lambda_bcs, omega_D, E_F)
        
        T_c = bcs.critical_temperature()
        print(f"\n临界温度: T_c = {T_c:.2f} K")
        
        # T = 0 时的性质
        Delta_0 = bcs.energy_gap(0)
        xi_0 = bcs.coherence_length(0)
        lambda_L_0 = bcs.penetration_depth(0)
        
        print(f"\nT = 0 时:")
        print(f"  能隙: Δ(0) = {Delta_0:.4e} eV")
        print(f"  相干长度: ξ(0) = {xi_0:.2e} m")
        print(f"  穿透深度: λ_L(0) = {lambda_L_0:.2e} m")
        
        # 温度依赖
        print(f"\n温度依赖:")
        temperatures = [0, 0.5*T_c, 0.9*T_c]
        for T in temperatures:
            Delta_T = bcs.energy_gap(T)
            H_c = bcs.critical_field(T)
            print(f"  T = {T:.2f} K:")
            print(f"    Δ(T) = {Delta_T:.4e} eV")
            print(f"    H_c(T) = {H_c:.2e} T")
        
        # 态密度
        print(f"\n超导态密度:")
        E_range = np.linspace(-3*Delta_0, 3*Delta_0, 100)
        dos_s = bcs.dos_superconductor(E_range, temperature=0)
        print(f"  能隙内 DOS: 0")
        print(f"  能隙边缘 DOS: 发散")
    
    @staticmethod
    def demonstrate_magnetism():
        """演示磁性"""
        print("\n" + "=" * 70)
        print("Heisenberg 自旋模型")
        print("=" * 70)
        
        # 铁磁耦合
        J_ferro = 1.0
        heisenberg = HeisenbergModel(J_ferro, lattice_size=100)
        
        print(f"\nHeisenberg 模型 (J={J_ferro}):")
        
        E = heisenberg.energy()
        M = heisenberg.magnetization()
        M_magnitude = np.linalg.norm(M)
        
        print(f"  能量: E = {E:.2f}")
        print(f"  磁化强度: |M| = {M_magnitude:.2f}")
        print(f"  是否铁磁: {heisenberg.is_ferromagnetic()}")
        
        # Curie 温度
        T_c = heisenberg.curie_temperature(S=0.5)
        print(f"  Curie 温度: T_c = {T_c:.2f}")
        
        # 反铁磁耦合
        print(f"\n反铁磁耦合 (J<0):")
        J_antiferro = -1.0
        print(f"  相邻自旋反平行")
        print(f"  Néel 温度: T_N ∝ |J|")
    
    @staticmethod
    def demonstrate_topology():
        """演示拓扑相"""
        print("\n" + "=" * 70)
        print("拓扑相与拓扑不变量")
        print("=" * 70)
        
        # 量子霍尔效应
        print(f"\n整数量子霍尔效应:")
        for nu in [1, 2, 3]:
            qhe = QuantumHallEffect(nu)
            sigma_xy = qhe.hall_conductivity()
            print(f"  ν = {nu}: σ_xy = {sigma_xy:.4e} S")
        
        print(f"\n分数量子霍尔效应:")
        for nu in [1/3, 2/5, 5/2]:
            qhe = QuantumHallEffect(nu)
            sigma_xy = qhe.hall_conductivity()
            print(f"  ν = {nu:.3f}: σ_xy = {sigma_xy:.4e} S")
        
        # 拓扑绝缘体
        print(f"\n拓扑绝缘体:")
        ti = TopologicalInsulator(bulk_gap=0.3, surface_states=True)
        print(f"  体能隙: {ti.E_g} eV")
        print(f"  表面态: {ti.has_surface_states}")
        print(f"  Z₂ 不变量: ν = {ti.z2_invariant()}")
        print(f"  拓扑性质: {'非平凡' if ti.is_topological() else '平凡'}")
        
        print(f"\n拓扑保护:")
        print(f"  边界态不受杂质散射破坏")
        print(f"  体现体-边对应原理")


def main():
    """主函数：完整演示"""
    
    print("=" * 70)
    print("DaoMath Condensed Matter Physics")
    print("从 ±0 构造凝聚态理论")
    print("=" * 70)
    
    # 1. 晶格结构
    CondensedMatterDemo.demonstrate_lattice()
    
    # 2. 能带理论
    CondensedMatterDemo.demonstrate_band_structure()
    
    # 3. 超导
    CondensedMatterDemo.demonstrate_superconductivity()
    
    # 4. 磁性
    CondensedMatterDemo.demonstrate_magnetism()
    
    # 5. 拓扑相
    CondensedMatterDemo.demonstrate_topology()
    
    print("\n" + "=" * 70)
    print("凝聚态物理构造完成！")
    print("=" * 70)


if __name__ == "__main__":
    main()
